home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / fortran / extenpr.zip / EXTEND.FOR < prev    next >
Text File  |  1988-06-04  |  16KB  |  709 lines
  1. C
  2. C    This is a collection of subroutines that perform extended precision
  3. C    floating point arithmetic.
  4. C
  5. C    By  Mark P. Esplin
  6. C        16 Shawsheen Rd.
  7. C        Billerica, Ma 01821    
  8. C
  9. C    Whole programs do not have to be written using extended
  10. C    precision.  Only the part that needs extra precision need use it,
  11. C    then convert back to regular REAL*8 numbers.
  12. C    The routines for addition and multiplication use the same algorithms 
  13. C    that are taught in grade school but using base 10000 numbers.
  14. C    This is not the most  efficient way to do it, but
  15. C    it makes it possible to use standard FORTRAN77 (almost).  
  16. C    These routines should work with any computer and compiler that supports
  17. C    INTEGER*2 and INTEGER*4 arithmetic.
  18. C    Each number is stored in an INTEGER*2 array of length N.
  19. C    The first word, IA(1), is the exponent in base 10000 digits.
  20. C    The least significant digit is in IA(2) and the most significant
  21. C    digit is in IA(N).  If IA(N) > 5000 the number is negative.
  22. C    When IA(N) = 5000 the number is assumed invalid.
  23. C    Base 10000 digits have worse roundoff characteristics than binary.
  24. C    Since roundoff causes the last base 10000 digit to be inaccurate, 
  25. C    that really means the last 4 decimal digits are inaccurate.  Be
  26. C    sure to use more digits than is necessary to allow for roundoff.
  27. C    N must be at least 6 or some routines won't work correctly.
  28. C
  29. C    Subroutines with their FORTRAN equivalent in this package are:
  30. C
  31. C    ADD(A,B,C,N)      C=A+B
  32. C    NEG(A,N)          A=-A
  33. C    MUL(A,B,C,N)      C=A*B
  34. C    INV(A,B,N,W)      B=1/A      W is work area of size 3*N
  35. C    EXSQRT(A,B,N,W)   B=SQRT(A)  W is work area of size 4*N
  36. C    EXHALF(A,N)       A=1/2      much faster than using INV
  37. C    PRTEXT(A,N)       WRITE(*,*) A
  38. C    WRTEXT(IU,A,N)    WRITE(IU,*) A
  39. C    CVI2EX(I2,A,N)    A=I2       I2 is INTEGER*2
  40. C    CVR8EX(R,A,N)     A=R        R is REAL*8
  41. C    CVEXR8(A,N,R)     R=A        R is REAL*8
  42. C    COPY(A,B,N)       B=A
  43. C    IEXCMP(A,B,N,W)    Compares two positive numbers A and B returns 1, 0, -1
  44. C    CHGPRE(A,NA,B,NB)  Changes the precision of A form NA to NB
  45. C
  46. C    There are other internal routines, but they would probably be of limited
  47. C    value for most users.
  48. C
  49. C
  50.     SUBROUTINE INCEX(IA,N)
  51.     INTEGER*2 IA(*),N
  52. C    This routine increments the mantissa of extended number IA.
  53. C    This routine is used by other routines to perform rounding.
  54.     INTEGER*2 IT,ICAR,I
  55.     ICAR=1
  56.     DO 10 I=2,N
  57.       IT=IA(I)+ICAR
  58.       ICAR=IT/10000
  59.       IA(I)=IT-ICAR*10000
  60. C    In almost all cases this routine will end after a few passes.
  61.     IF(ICAR.EQ.0) GOTO 20
  62. 10    CONTINUE
  63. C    Shift back if necessary
  64. C    The only time that a shift is necessary is when IA(N) = 5000 and
  65. C    the rest of the mantissa is zero.
  66. 20    CONTINUE
  67.     IF(IA(N).EQ.5000) THEN
  68.       IA(1)=IA(1)+1
  69.       DO 30 I=2,N-1
  70.         IA(I)=IA(I+1)
  71. 30      CONTINUE
  72.       IA(N)=0
  73. C      Don't need to round again since if a shift was necessary
  74. C      the least significant digit is zero.
  75.     ENDIF
  76.     RETURN
  77.     END
  78. C
  79. C
  80. C
  81. C
  82. C
  83.     FUNCTION ICKCON(IA,N)
  84. C    This routine is used to check if an extended number IA is converging
  85. C    to an integer value. 
  86. C    This routine would probably not be very useful for the average user.
  87. C    The value returned is 1 if IA has converged.  The list significant digit
  88. C    is not checked since it might be just round off.
  89. C    This function return a 1 in IA has converged to an integer otherwise
  90. C    it returns a 0.
  91.     INTEGER*2 IA(*)
  92.     INTEGER*2 N,I
  93.     INTEGER*2 ITEST
  94.     ICKCON=0
  95. C    The number must have already converged past the N-2 th digit before
  96. C    this routine is called.
  97.     ITEST=IA(N-2)
  98.     DO 10 I=3,N-3
  99.       IF(IA(I).NE.ITEST) RETURN
  100. 10    CONTINUE
  101.     ICKCON=1
  102.     RETURN
  103.     END
  104.     SUBROUTINE NEG(IA,N)
  105. C    This routine changes the sign of extended precision number IA.
  106.     INTEGER*2 IA(*),N
  107.     INTEGER*2 I,J
  108.     J=2
  109.     DO 10 I=2,N
  110.       IF(IA(J).NE.0) GOTO 20
  111.       J=J+1
  112. 10    CONTINUE
  113. C    Number was zero
  114.     RETURN
  115. C    The first non-zero digit has been reached
  116. 20    IA(J)=10000-IA(J)
  117.     J=J+1
  118.     DO 30 I=J,N
  119.       IA(I)=9999-IA(I)
  120. 30    CONTINUE
  121.     RETURN
  122.     END
  123. C
  124. C
  125. C
  126. C
  127.     SUBROUTINE PRTEXT(IA,N)
  128. C    This routine writes extended precision numbers to output.
  129. C    This routine could use some work to make a more readable output.
  130.     INTEGER*2 IA(*),N
  131.     INTEGER*2 I,INEG
  132.     INEG=0
  133.     WRITE(*,*) ' '
  134.     IF(IA(N).GT.5000) THEN
  135.       WRITE(*,*) 'Minus'
  136.       CALL NEG(IA,N)
  137.       INEG=1
  138.     ENDIF
  139.     WRITE(*,900) IA(1)
  140. 900    FORMAT(' Exponent  ',I6)
  141.     WRITE(*,910) (IA(I),I=N,2,-1)
  142. 910    FORMAT(10I5)
  143. C    Restore IA to previous value.
  144.     IF(INEG.NE.0) CALL NEG(IA,N)
  145.     RETURN
  146.     END
  147. C
  148. C
  149. C
  150. C
  151.     SUBROUTINE WRTEXT(IUNIT,IA,N)
  152. C    This routine writes extended precision numbers to file IUNIT.
  153. C    This routine could use some work to make a more readable output.
  154.     INTEGER*2 IUNIT,IA(*),N
  155.     INTEGER*2 I,INEG
  156.     INEG=0
  157.     WRITE(IUNIT,*) ' '
  158.     IF(IA(N).GT.5000) THEN
  159.       WRITE(IUNIT,*) 'Minus'
  160.       CALL NEG(IA,N)
  161.       INEG=1
  162.     ENDIF
  163.     WRITE(IUNIT,900) IA(1)
  164. 900    FORMAT(' Exponent  ',I6)
  165.     WRITE(IUNIT,910) (IA(I),I=N,2,-1)
  166. 910    FORMAT(10I5)
  167. C    Restore IA to previous value.
  168.     IF(INEG.NE.0) CALL NEG(IA,N)
  169.     RETURN
  170.     END
  171. C
  172. C
  173. C
  174. C
  175. C
  176.     SUBROUTINE NORMIZ(IA,N)
  177. C    This routine is used by other routines to normalize
  178. C    floating point number after addition and other such processes.
  179.     INTEGER*2 IA(*),N
  180.     INTEGER*2 I,ISHFT
  181.     IF(IA(N).GT.5000) THEN
  182. C      This is a negative number.
  183.       DO 10 I=N,2,-1
  184.         IF(IA(I).NE.9999) GOTO 20
  185. 10      CONTINUE
  186.       I=I+1
  187. 20      ISHFT=N-I
  188.       IF(IA(I).LE.5000) ISHFT=ISHFT-1
  189.     ELSE
  190. C      This is a positive number.
  191.       DO 30 I=N,2,-1
  192.         IF(IA(I).NE.0) GOTO 40
  193. 30      CONTINUE
  194. C      This number is zero
  195.       RETURN
  196. 40      ISHFT=N-I
  197.       IF(IA(I).GE.5000) ISHFT=ISHFT-1
  198.     ENDIF
  199.     IF(ISHFT.EQ.0) RETURN
  200.     IA(1)=IA(1)-ISHFT
  201. C    Shift left the appropriate amount.
  202.     DO 50 I=N,ISHFT+2,-1
  203.       IA(I)=IA(I-ISHFT)
  204. 50    CONTINUE
  205. C    Zero unknown digits.
  206.     DO 60 I=2,ISHFT+1
  207.       IA(I)=0
  208. 60    CONTINUE
  209.     RETURN
  210.     END
  211. C
  212. C
  213. C
  214. C
  215. C
  216.     SUBROUTINE ADD(IA,IB,IC,N)
  217. C    This routine performs extended precision addition.
  218. C    IC = IA + IB. IA and IB are not changed.
  219.     INTEGER*2 IA(*),IB(*),IC(*)
  220.     INTEGER*2 N,IPA,IPB,NML,I
  221.     INTEGER*2 ISGNA,ISGNB,ISGNC,ISGEXT
  222.     INTEGER*2 IT,ICAR
  223.     ISGNA=1
  224.     IF(IA(N).GT.5000) ISGNA=-1
  225.     ISGNB=1
  226.     IF(IB(N).GT.5000) ISGNB=-1
  227. C    Figure out pointers for main loop
  228.     IF(IA(1).GE.IB(1)) THEN
  229.       IC(1)=IA(1)
  230.       IPA=2
  231.       IPB=IA(1)-IB(1)+2
  232.       NML=N-IPB+2
  233.       IF(NML.LT.2) THEN
  234. C        The difference between the two numbers is too large.
  235.         CALL COPY(IA,IC,N)
  236.         RETURN
  237.       ENDIF
  238.       IF(IPB.EQ.2) THEN
  239.         ICAR=0
  240.       ELSE
  241.         ICAR=(IB(IPB-1)+5000)/10000
  242.       ENDIF
  243.     ELSE
  244.       IC(1)=IB(1)
  245.       IPB=2
  246.       IPA=IB(1)-IA(1)+2
  247.       NML=N-IPA+2
  248.       IF(NML.LT.2) THEN
  249. C        The difference between the two numbers is too large.
  250.         CALL COPY(IB,IC,N)
  251.         RETURN
  252.       ENDIF
  253.       ICAR=(IA(IPA-1)+5000)/10000
  254.     ENDIF
  255. C    Main addition loop.
  256.     DO 10 I=2,NML
  257.       IT=IA(IPA)
  258.       IT=IT+IB(IPB)+ICAR
  259.       ICAR=IT/10000
  260.       IC(I)=IT-ICAR*10000
  261.       IPA=IPA+1
  262.       IPB=IPB+1
  263. 10    CONTINUE    
  264. C    Do sign extended part.
  265.     IF(IA(1).GT.IB(1)) THEN
  266.       IF(ISGNB.GE.1) THEN
  267.         ISGEXT=0
  268.       ELSE
  269.         ISGEXT=9999
  270.       ENDIF
  271.       DO 20 I=NML+1,N
  272.         IT=IA(I)
  273.         IT=IT+ISGEXT+ICAR
  274.         ICAR=IT/10000
  275.         IC(I)=IT-ICAR*10000
  276. 20      CONTINUE
  277.     ELSE
  278.       IF(ISGNA.GE.1) THEN
  279.         ISGEXT=0
  280.       ELSE
  281.         ISGEXT=9999
  282.       ENDIF
  283.       DO 30 I=NML+1,N
  284.         IT=IB(I)
  285.         IT=IT+ISGEXT+ICAR
  286.         ICAR=IT/10000
  287.         IC(I)=IT-ICAR*10000
  288. 30      CONTINUE
  289.     ENDIF  
  290. C    Shift right one place if it has overflowed during addition.
  291.     IF(ISGNA*ISGNB.EQ.1) THEN
  292.       ISGNC=1
  293.       IF(IC(N).GE.5000) ISGNC=-1
  294.       IF(ISGNA*ISGNC.EQ.-1) THEN
  295. C        An overflow has occurred.
  296.         IC(1)=IC(1)+1
  297.         ICAR=IC(2)
  298.         DO 40 I=2,N-1
  299.           IC(I)=IC(I+1)
  300. 40        CONTINUE
  301.         IC(N)=ISGEXT
  302. C        Do rounding and return.
  303.         IF(ICAR.GE.5000) CALL INCEX(IC,N)
  304.         RETURN
  305.       ENDIF
  306.     ENDIF
  307. C    Normalize answer.
  308.     CALL NORMIZ(IC,N)
  309.     RETURN
  310.     END
  311. C
  312. C
  313. C
  314. C
  315.     SUBROUTINE MUL(IA,IB,IC,N)
  316. C    This routine multiplies two extended precision numbers.
  317. C    N must be at least 4 or this routine won't work.
  318. C    IC = IA*IB
  319.     INTEGER*2 IA(*),IB(*),IC(*),N
  320.     INTEGER*2 ISNA,ISNB,IXC1,IXC2,ISHFT
  321. C    The computation is carried out with two extra placesIXC1 and IXC2
  322. C    This makes it so precision is not lost if the